Fast minimum variance wavefront reconstruction for extremely large telescopes 
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We present a new algorithm, FRiM (FRactal Iterative Method), aiming at the recon- 
struction of the optical wavefront from measurements provided by a wavefront sensor. 
As our application is adaptive optics on extremely large telescopes, our algorithm was 
designed with speed and best quality in mind. The latter is achieved thanks to a reg- 
ularization which enforces prior statistics. To solve the regularized problem, we use 
the conjugate gradient method which takes advantage of the sparsity of the wavefront 
sensor model matrix and avoids the storage and inversion of a huge matrix. The prior 
covariance matrix is however non-sparse and we derive a fractal approximation to the 
Karhunen-Loeve basis thanks to which the regularization by Kolmogorov statistics 
can be computed in 0{N) operations, being the number of phase samples to es- 
timate. Finally, we propose an effective preconditioning which also scales as 0{N) 
and yields the solution in 5-10 conjugate gradient iterations for any A^. The result- 
ing algorithm is therefore 0{N). As an example, for a 128 x 128 Shack-Hartmann 
wavefront sensor, FRiM appears to be more than 100 times faster than the classical 
vector-matrix multiplication method. 
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1. Introduction 

The standard and most used method for adaptive optics (AO) 
control is based on a vector-matrix multiply (VMM) of the 
vector of wavefront sensor measurements by the so-called 
control matrix [1]. This operation gives an update of the com- 
mands to be sent to the deformable mirrors to adjust the cor- 
rection of the corrugated incoming wavefronts. The control 
matrix is precomputed, generally using modal control opti- 
mization [2]. The complexity of computing the control ma- 
trix using standard methods scales as 0{N^), where is the 
number of unknowns (phase samples or actuator commands), 
and applying real time VMM scales as 0{N~). This compu- 
tational burden can be reasonably handled on current AO sys- 
tems where < lO-'. 

For future Extremely Large Telescopes (ELT's), the number 
of actuators beeing considered is in the range lO"* - 10^. This 
huge increase is the result of both the larger diameter of the 
ELTs [3] and the emergence of new architectures for the AO 
systems, using either a greater density of actuators (Extreme 
AO) or combining several deformable mirrors and wavefront 
sensors (multi-conjugate AO, multi-object AO) [4]. The nec- 
essary computational power for real time control on such sys- 
tems is currently unattainable when using standard methods. 

More efficient algorithms are thus required and have been 
developed in recent years. Poyneer et al. [5] have derived an 
accurate Fourier transform wavefront reconstructor by solv- 
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ing the boundary problem in circular apertures. This recon- 
structor scales as 0{N log A^) and is shown to be effective for 
Extreme AO [6]. MacMartin [7] studied several approximate 
approaches such as a multiple-layer hierarchic reconstruction, 
which scales as 0{N). 

Although least-squares algorithms give suitable results for 
single star AO systems (classical on-axis AO or Extreme AO), 
minimum variance reconstruction is required to minimize the 
eff'ects of the missing data or unseen modes in the other AO 
schemes [8]. In the context of minimum variance for multi- 
conjugate AO, Ellerbroek [9] could apply sparse matrix tech- 
niques (Cholesky factorization) using a sparse approximation 
of the turbulence statistics, and introducing as low-rank ad- 
justments, the nonsparse matrix terms arising from the global 
tip/tilt measurement errors associated with laser guide stars. 
However the interactions between the layers in their tomo- 
graphic modeling reduce the efficiency of the sparse direct 
decomposition methods [10]. 

Iterative methods are also extensively studied in this con- 
text. Their main asset is their ability to iteratively compute 
the unknowns from the measurements using direct sparse ma- 
trices, and so the storage of a precomputed inverse full matrix 
is not necessary. One major problem with iterative methods 
is the increase in the number of iterations with the number of 
unknowns to estimate [11-13]. As an example. Wild et al. 
[ 14] have proposed to use the closed-loop AO system itself as 
an iterative processor, but the performance of the least squares 
reconstruction depends on the loop frequency of the AO sys- 
tem, which should be higher than usual. 

The most successful iterative methods in AO are now based 
on preconditioned conjugate gradients (PCG) [15], where 
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some of the previous approximate reconstruction methods are 
embedded as preconditioners to ensure a small number of it- 
erations (see section 4). Gilles et al. [16] have described a 
multigrid PCG algorithm, mainly aimed at Extreme AO and 
scaling as 0{N log N). The multigrid preconditioner is some- 
what related to the multiple-layers hierarchic reconstruction 
[7]. This wavefront reconstruction method has been improved 
with a faster approximation to the turbulence statistics, scal- 
ing as 0{N) [17]. The multigrid PCG algorithm has also been 
developed for multi-conjugate AO [18]. In this case, the struc- 
ture of the matrix is more complex and brings some limita- 
tions. More recently, a Fourier domain preconditioner was in- 
troduced [19, 20] in the context of multi-conjugate AO, with 
a faster reconstruction than multigrid PCG. In this case, the 
preconditioner is related to the Fourier transform wavefront 
reconstructor [5]. Both multigrid and Fourier domain pre- 
conditioners were examined for the Thirty Meter Telescope 
project [21, 22]. 

In this work, we propose novel methods to address the two 
critical points previously seen in iterative methods for wave- 
front reconstruction: estimation of the atmospheric phase co- 
variance matrix and preconditioning. 

We need a sparse representation of the inverse of the atmo- 
spheric phase covariance matrix to efficiently introduce pri- 
ors in the minimum variance estimator. Currently, we can 
choose between a good representation in the Fourier domain 
with 0{N log A^) complexity [16, 19] and a widely used sparse 
biharmonic approximation introduced by Ellerbroek [9], less 
accurate [19], but scaUng as 0{N). With FRiM, we introduce 
a so-called "fractal operator" as a multiscale algorithm with 
0{N) complexity. This operator, both accurate and very fast, 
was inspired by the mid-point method of Lane et al. [23] to 
generate a Kolmogorov phase screen. It can be used for any 
wavefront structure function. It allows us to very efficiently 
apply the inverse of the phase covariance matrix to any vector 

We show that this fractal operator is also very efficient when 
used as a preconditioner. It allows the wavefront reconstruc- 
tion to be iteratively computed in a space of statistically inde- 
pendent modes. We additionally use a classical Jacobi precon- 
ditioner, or a new "optimal diagonal preconditioner" to further 
improve the convergence. The number of iterations is < 10 for 
a fuU wavefront reconstruction whatever the size of the sys- 
tem, with a number of floating point operations ~ 34 x per 
iteration. The method is therefore globally 0{N). 

In the following, we first derive the analytical expression 
for the minimum variance restored wavefront and the equa- 
tions to be solved. We then introduce the fractal operator al- 
lowing fast computation of the regularization term in an itera- 
tive method such as conjugate gradients. We then propose two 
fast preconditioners to further speed up the iterative algorithm. 
We finally use numerical simulations to test the performances 
of FRiM. 

2. Minimum variance solution 

A. Model of data 

We assume that the wavefront sensor provides measurements 
of spatial derivatives (slopes or curvatures) of the phase. 



which are linearly related to the wavefront seen by the sen- 
sor: 

d^Sw + n (1) 

where d e R*' is the data vector provided by the sensor, w e 
R'^ is the vector of sampled wavefront values, S e R*^^'^ is the 
sensor response matrix and n e R*' accounts for the noise and 
model errors. This equation is general as long as the wavefront 
sensor is linear. As a typical case, we will however consider a 
Shack-Hartmann wavefront sensor with Fried geometry [24] 
in our simulations and for the evaluation of the efficiency of 
the algorithms. 

B. Optimal wavefront reconstructor 

The estimation of the wavefront w given the data d is an in- 
verse problem which must be solved using proper regular- 
ization in order to improve the quality of the solution while 
avoiding noise amplification or ambiguities due to missing 
data [25]. In order to keep the problem as simple as possi- 
ble, we first introduce the requirement that the solution be a 
linear function of the data, i.e. the restored wavefront satisfies: 



where R is the restoration matrix and d the wavefront sensor 
measurements. Some quality criterion is needed to derive the 
expression for the restoration matrix R. For instance, we can 
require that, on average, the difference between the restored 
wavefront w and the true wavefront w be as small as possi- 
ble by minimizing (||h' - w\\~} where (■) denotes the expected 
value of its argument. It is interesting to note that minimizing 
(on average) the variance of the residual wavefront yields the 
optimal Strehl ratio [26] since: 

SR-exp(-4 r [w{r)-w{r)f dr] (3) 

\ Jpupil / 

where r is the position in the pupil, ^ is the area of the pupil 
and w (r) is the wavefront phase in radian units. The best re- 
construction matrix according to our criterion then satisfies: 

R"'' ^argmindlR-rf-vfll"). (4) 

R 

Accounting for the facts that the wavefront w and the errors 
n are uncorrected and have zero means, i.e. («) = and 
{w} = 0, the minimum variance reconstructor expands as [27]: 

R"'" = C ■ ■ (S ■ C„, ■ C„)"' , (5) 

where C„ {n ■ n^} is the covariance matrix of the errors and 

Cw = {w ■ w^) is the a priori covariance matrix of the wave- 
front. Applying this reconstructor to the data d requires solv- 
ing a linear problem with as many equations as there are mea- 
surements. Generally, wavefront sensors provide more mea- 
surements than wavefront samples (about twice as many for a 
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Shack-Hartmann or a curvature sensor). Fortunately, from the 
following obvious identities [28]: 

S ■ C|j ' S ■ C^; ■ S ~t" S — S ■ C|| ' ■ Cw ' S + 

= (sT.c„'-s + C)-c,-sT, 

we can rewrite the optimal reconstructor in Eq. (5) as: 

R"'" = (ST-C„'-S + C,'f'-ST-C„' (6) 

which involves solving just as many linear equations as there 
are wavefront samples. The linear reconstructor defined in 
Eq. (6) is the expression to be preferred in our case. 

C. Links with other approaches 

Using Eq. (6) for the reconstructor, the minimum variance re- 
stored wavefront is given by: 

R"'' ■ = (S^ ■ C„' ■ S + C^')"' ■ S'^ ■ C„' ■ 
which is also the solution of the quadratic problem: 
H-"'' = argmin{(S ■ w - df ■ ■ {S ■ w - d) + ■ ■ h-} 

w 

where (S • w - d)^ ■ C„^ ■ (S ■ w - rf) is the so-called which 
measures the discrepancy between the data and their model 
and w"^ ■ C^' • w is a Tikhonov regularization term which en- 
forces a priori covariance of the unknowns. Thus Eq. (6) is 
also the result of the maximum a posteriori (MAP) problem. 
Here, the usual hyper-parameter is hidden in C„, which is pro- 
portional to (Dlrof^^, where ro is the Fried parameter [29]. 
As already noted by other authors (see e.g. Rousset [30]), the 
minimum variance estimator is directly related to Wiener op- 
timal filtering. 

Actual adaptive optics systems make use of some expansion 
of the wavefront on a basis of modes, regularization being 
achieved by setting the ill-conditioned modes to zero. This 
technique is similar to truncated singular value decomposition 
(TSVD) [30]. Since truncation results in aliasing, we expect 
that the MAP solution will be a better approximation to the 
wavefront. 

D. Iterative Method 

The optimal wavefront can be computed in different ways. For 
instance, the matrix R can be computed once, using Eq. (5) or 
Eq. (6), and then applied to every data set d. Since it requires 
the numerical inversion of an NxN matrix, the direct compu- 
tation of R scales as 0{N^) operations [13]. The reconstructor 
R is a X M matrix and is not sparse in practice. Hence, the 
storage of R requires M N ~ 2N^ floating point numbers and 
computing R ■ d requires a: 2M N x 4N^ floating point oper- 
ations. For large numbers of degrees of freedom oc (D/ro)^, 
the computer time spent by the matrix-vector multiplication 
can be too long for real time applications. Moreover the mem- 
ory requirement (e.g. for ^ 10"*, 1.5 Gb of memory are 
needed to store R) may be such that memory page faults dom- 
inate the computation time of matrix-vector multiplication. 



initialisation: 

compute rg = b - A ■ Xo for some initial guess Xq 

let = 
until convergence do 

solve M • Zj. = Tj. for Zk (apply preconditioner) 

Pk = rj -Zk 

ifyt = 0, then 

Pk = 
else 

Pk = Zk + (PklPk-l)Pk-l 
endif 

9a = A • Pk 

ak = Pk/ipJ-qk) (optimal step size) 

Xk+l = Xk + CfkPk 

fk+i =rk-ak Qk 
done 



Fig. 1. Preconditioned conjugate gradient algorithm for solv- 
ing A-JC = b where A is a symmetric positive definite ma- 
trix and M is a preconditioner The unpreconditioned version 
of the algorithm is simply obtained by taking M = I, hence 
Zk = rk. 

In order to avoid the direct matrix inversion and the matrix- 
vector product required by the explicit computation of R, we 
use an iterative method to solve the linear system 

(sT.C„'-S + C)-H' = ST-C„'-rf (7) 

which leads to the optimal wavefront w for every data set d. 
For the purpose of the discussion, Eq. (7) can be put in a more 
generic form: 

A x^b (8) 

where, in the case of Eq. (7), x - w and: 

A = ■ C„' ■ S + C,' (9) 

is the so-called left hand side matrix, whereas 

b^S^C„^d (10) 

is the so-called right hand side vector 

Barrett et al. [15] have reviewed a number of iterative al- 
gorithms for solving linear systems like (8). An advantage of 
these methods is that they do not explicitly require the ma- 
trix A; it is sufficient to be able to compute the product of 
matrix A (or its transpose) with any given vector The itera- 
tive algorithm therefore fully benefits from the possibility to 
compute the matrix-vector products in much less than 0{N~) 
operations when A is sparse or has some special structure. 
This is particularly relevant in our case since applying A can 
be achieved by matrix-vector products very fast to compute 
as shown in Sect. 2E and Sect. 2F. The drawback of iter- 
ative methods is that the computational burden scales as the 
number of iterations required to approximate the solution with 
sufficient precision. In the worst case, the number of itera- 
tions can theoretically be as high as the number of unknowns 
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w{x, y+ a) 




w{x + a, y + a) 



dy{x, y) 



dxU y) 



w(x" + a, y) 



Fig. 2. Wavefront sensor with Fried geometry as used for our 
simulations. The black circles stand for phase samples w{x, y), 
at the corners of the square subapertures of size a. This model 
is exact if we assume that the wavefront at any point in the 
pupil is obtained from a bilinear interpolation of phase sam- 
ples at the corner of the subapertures. 



[12, 13]. In practice and because of numerical rounding 
errors, ill-conditioning of the system in Eq. (7) can result in 
a much higher number of iterations, even on small systems. 
This problem can however be greatly reduced by means of a 
good preconditioner [ 12, 15]. 

By construction, A given by Eq. (9) is a symmetric posi- 
tive definite matrix and the conjugate gradient (CG) [15] is 
the iterative method of choice to solve the system in Eq. (8). 
Figure 1 shows the steps of the CG algorithm to solve the sys- 
tem A x = b. This method is known to have a super-linear 
rate of convergence [12], and can be accelerated by using a 
proper preconditioner M a; A for which solving M ■ z = r for 
z (with r - b - A ■ x) is much cheaper than solving Eq. (8) 
for X. The preconditioner can also be directly specified by its 
inverse Q = M"' such that Q a; A"' and then z = Q ■ r in the 
CG algorithm. Without a preconditioner, taking M = Q = I, 
where I is the identity matrix, yields the unpreconditioned ver- 
sion of the CG algorithm. In Sect. 4 we investigate various 
means to obtain an eff'ective preconditioner for the wavefront 
reconstruction problem. 

In the remainder of this section, we derive means to quickly 
compute the dot product with the matrix A in Eq. (9). To that 
end, we consider separately the Hessian matrix ■ C„' ■ S of 
the likelihood term and that of the regularization term C^' . 

E. Computation of the likelihood term 

Most adaptive optics systems use either a Shack-Hartmann 
sensor which provides measurements of the local gradient of 
the wavefront or a curvature sensor which measures the local 
curvature of the wavefront [1]. Since such sensors probe lo- 
cal spatial derivatives of the wavefront, their response can be 
approximated by local finite differences which yields a very 
sparse linear operator S. Though some non-sparse matrix 
terms can appear due to tilt indetermination with laser guide 
stars or to take account of natural guide star tip/tilt sensors. 
Owing to the low rank of these modes, sparse matrix models 
can still be applied [9]. Thus, denoting A^^if the number of 
wavefront samples required to compute the local finite differ- 
ences, only a M X A^dif out of M x N coefficients of S are 



non-zero. For instance. Fig. 2 shows the Fried geometry of 
the Shack-Hartmann sensor model [24] which we used in our 
numerical simulations. The error free slopes are related to the 
wavefront by: 



dx{x,y) - J [w{x + a,y + a) + w(x + a,y) 
- w(x, y + a) - w(x, y)] 

dy{x,y) — J [w{x + a,y + a) — w(x + a,y) 
+ w(x, y + a) - w(x, y)] 



(11) 



where (x, y) are the pupil coordinates, d^ and dy are the slopes 
along the x and y directions and a is the sampling step. Hence 
Ndif - 4, in our case, whatever the number of degrees of free- 
dom. Besides, to a good approximation, wavefront sensors 
provide uncorrected measurements [1], hence the covariance 
matrix C„ of the errors can be taken as a diagonal matrix: 



C„ a diag(Var(n)) 



(12) 



where Var(M) is the vector of noise and error variances. Since 
C„ is diagonal, its inverse C„' is diagonal and trivial to com- 
pute. Finally, the matrices S and C„' are sparse and the dot 
product by ■ C„' ■ S can be therefore computed in 0{N) 
operations. 

F. Fast estimation of the regularization term 

Unlike C„ and C„', neither C», nor C^' is sparse. We intro- 
duce here a way to derive an approximation for C^' so that it 
can be applied to a vector with a small number of operations. 
We first consider the following decomposition of C,^: 



Ch- — K ■ 



(13) 



where K is a square invertible matrix. Since is positive 
definite, there exists a number of possibilities for such a fac- 
torization: Cholesky decomposition [9, 13, 17], spectral fac- 
torization, etc. We then use this decomposition to define new 
variables u based on the wavefront w: 



def 1 



(14) 



The expected value of u is: (m) 
ance matrix therefore satisfies: 



K ' ■ (h") = and its covari- 



C„ = {u ■ u^) = ■ {w ■ w^) ■ 
= K ' C„, K^ = I, 

which shows that the new variables are independent and iden- 
tically distributed following a normal law: u ~ N{0, 1). This 
gives rise to a method for generating wavefronts since from 
a set u of N independent random values following a normal 
law, taking w = K - u yields a random wavefront with the ex- 
pected covariance. Finally, using this re-parametrization, it is 
possible to rewrite the regularization term as: 



(15) 



Then, depending on whether the problem is solved for the 
wavefront samples w or for the so-called wavefront genera- 
tors u (cf equations (38) and (39) in Sect. 5), each conjugate 



gradient iteration would be cheap to compute providing either 
(i) operators K"' and K"^ are fast to apply, or (ii) operators K 
and are fast to apply. 

A comparable re-parametrization has been introduced by 
Roddier [31] for generating turbulent wavefronts using a 
Zernike expansion of randomly weighted Karhunen-Loeve 
functions. This however requires to diagonalize a huge N xN 
matrix C^., a procedure that costs at least 0{N^) operations, 
and would give a slow operator K (or K"') taking 0{N^) op- 
erations to apply. 

Exploiting the fractal structure of turbulent wavefronts. 
Lane et al. [23] have derived a fast method to generate wave- 
fronts by a mid-point algorithm. In what follows, we show 
that their method amounts to approximating the effect of op- 
erator K in 0{N) operations and we derive algorithms to apply 
the corresponding K K^, and K"^ operators that also take 
0{N) operations. We propose to use these so-called fractal op- 
erators for fast computation of the regularization and also as 
effective pre-conditioners to speed-up the conjugate-gradient 
iterations. 

3. Fractal operators 

A. Principle and structure function 

The mid-point algorithm [23] starts at the largest scales of 
the wavefront and step-by-step builds smaller scales by in- 
terpolating the wavefront values at the previous scale and by 
adding a random value with a standard deviation computed 
so that the new wavefront values and their neighbors have the 
expected structure function. Using to denote the linear op- 
erator which generates the wavefront values at the y-th scale, 
the linear operator K can be factorized as: 



K — Ki ■ 



K„ 



(16) 



where p is the number of scales, generates the 4 outermost 
wavefront values and Ki generates the wavefront values at 
the finest scale. The original mid-point algorithm cannot be 
used directly for our needs because it is not invertible. In this 
section, we reconsider the mid-point algorithm to derive new 
expressions for the K/s such that they are sparse, invertible 
and such that their inverses are also sparse. 

The structure function of the wavefront is the expected 
value of the quadratic difference between two phases of a tur- 
bulent wavefront: 



<[w(r,)-w(r,)]'>=/(|r,-r,|). 



where, e.g.: 



f{r) = 6.88 X (r/ro)^ 



5/3 



(17) 



(18) 



for a turbulent wavefront obeying Kolmogorov's law. The 
structure function is stationary (shift-invariant) and isotropic 
since it only depends on the distance |r, - rj\ between the con- 
sidered positions r,- and rj in the wavefront. From the struc- 
ture function, we can deduce the covariance of the wavefront 
between two positions in the pupil: 




Fig. 3. The four initial values for wavefront generation, at the 
corners of the support. 



with Wi = w(r,) the wavefront phase at position r,, cr^ - 
Var(w,), and fj - /(|r; - rj\) the structure function between 
wavefront samples ; and j. The wavefront variances (thus the 
covariance) are not defined for pure Kolmogorov statistics but 
can be defined by other models of the turbulence such as the 
von Karman model. Nevertheless, any structure function / 
can be used by our algorithm: in case the variance is unde- 
fined, we will show that the crr's appear as free parameters 
and that choosing suitable variance values is not a problem. 

B. Generation of outermost values 

The first point to address is the initialization of the mid-point 
recursion, that is the generation of the four outermost corner 
values. Lane et al. [23] used 6 random values to generate the 4 
initial corners. It is however required to use exactly the same 
number of random values u as there are wavefront samples 
in w otherwise the corresponding linear operator K cannot be 
invertible. This is possible by slightly modifying their original 
algorithm. 

The four initial wavefront values (Fig. 3) have the following 
covariance matrix: 



Cout 



' C\ C2 Ci 

Ci Co Ci C2 

C2 Ci Co Ci 

( ci C2 ci Co ; 



with 



Co = cr- 

ci = (T^-f(D)l2 
C2 = o-2-/(V2D)/2 



where en is the variance (assumed to be the same) of the four 
initial phases and where D is the distance between points 1 
and 2 (see Fig. 3). Having the same variances en for the four 
outermost wavefront phases seems natural since none of these 
points plays a particular role. For the four outer wavefront 
samples, the matrix of eigenvectors of Cout is: 



Zniir — 



' 1/2 -1/2 

1/2 1/2 -I/V2 

1/2 -1/2 

1/2 1/2 I/V2 





-I/V2 




Note that the eigenvectors (columns) defined on these four 
samples are (in order) piston, waffle [7], tip and tilt. The eigen- 
values are: 



/Inut — 



(19) 



' Co + 2 Ci H- C2 
Co - 2 Ci H- C2 

Co - C2 
Co - C2 



4t^2-/(D)-/(V2D)/2^ 
/(D)-/(V2D)/2 
/(V2D)/2 
/(V2D)/2 



In the case of pure Kolmogorov statistics, cr^ must be cho- 
sen so that K is invertible. This is achieved if the eigenvalue 
of the piston-\\\ie mode is strictly positive, hence: 

(r2>/(D)/4 + /(V2D)/8. 

We have chosen so that the smallest covariance, which is 
c(V2Z)) between the most remote points, is exactly zero: 



(20) 



Of course, when a von Karman model of turbulence is chosen, 
both cr- and / are fixed by the model; Eq. (20) is to be used 
only for the Kolmogorov case. 

A possible expression for the operator Kom, such that Cout = 

Kout ■ 



K^m, is: 



Knnt — — 



a -b -c 0^ 
a b —c 
a —b c 
a b c 



with: a = ^J4c^^ - f(D) - f(^D)/2 , 
b= ylf{D)-f{^.D)l2, 

C= yffi^, 



(21) 




Fig. 4. Wavefront refinement. To generate a grid with cell size 
r/2, new values (in gray) are generated from wavefront values 
(in white) of a grid with cell size equal to r. Top left: new 
value from 4 values r/ V2 apart. Top right: new edge value 
from 3 values r/2 apart. Bottom: new value from 4 values r/2 
apart. 



from which K^^^j is: 



I /a 


1/fl 




1/fl 


l/b 


l/b 


-l/b 


l/b 


2/c 





2/c 








-2/c 





2/c 



(22) 



The operator Kp in Eq. (16) is obtained simply from Kouf 
Indeed, K,, is essentially the identity matrix except for 16 non- 
zero coefficients corresponding to the outermost corners and 
which are given by Kout- The same rules yield K^' from K^^'j. 

C. Generation of wavefront samples at smaller scales 

Given the wavefront with a sampling step r, the mid-point 
algorithm generates a refined wavefront with a sampling of 
r/2 using a perturbed interpolation: 



Wo = ao Mo ■ 



(23) 



where wq is the wavefront value at the mid-point position. 
Mo ~ A/'(0, 1) is a normally distributed random value and 
Nint is the number of wavefront samples from the previous 
scale which are used to generate the new sample (see Fig. 4). 
Equation (23) comes from a generalization of the principle of 
the original mid-point algorithm. Since we proceed from the 
largest scale to smaller ones, all the operations can be done in- 
place: the value of wq computed according to Eq. (23) replac- 
ing that of mq. In other words, the input and output vectors, u 
and w, can share the same area of the computer memory. It is 
then immediately apparent that a random wavefront computed 



by this algorithm scales as 0{N\at xN) - 0{N) since the num- 
ber of neighbors Nint ~ 4 does not depend on the number of 
wavefront samples A^. 

The A^int + 1 scalars aj have to be adjusted so that the struc- 
ture function between wq and any of the Wi^i^,,,^^.^^ matches the 
turbulence statistics: 

f.O <(W0 - Wif) 

Mm 



j=l l<j<k<Ni^, 



K k=l ) 



./=I 



(24) 



Note that, to obtain this equation, we have accounted for the 
fact that since mo ~ N {0, 1) and Wj=i,...,Arj„, are uncorrected, 
then (mq) - 1 and (mq Wj=i,...,iv.„,) - 0. The system (24) gives 
A^int equations, whereas there are Nim + 1 unknown parameters 
{q-o, . . . , o'A'ini): an additional constraint is needed. 

In the original mid-point algorithm. Lane et al. [23] choose 
to normalize the sum of the interpolation coefficients and use 
the constraint that ctj - 1 ■ In that case, Eq. (24) simpli- 
fies and the coefficients are obtained by solving: 



j=l l<j<k<Ni«, 

s.t. YjOj^l 



(25) 



Note that all the variances crj are implicit with this constraint. 
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We consider here another constraint which is to have the 
same variance, say cr^, for all the wavefront samples. In 
other words, we consider a wavefront with stationary (shift- 
invariant) statistical properties. This is justified by our objec- 
tive to reconstruct phase corrugations in several layers for at- 
mospheric tomography. Indeed, since the beams coming from 
different directions in the field of view are not superimposed 
in the layers, this condition allows the wavefront statistics to 
remain the same for all the beams. With this choice, the addi- 
tional equation is provided by (w^) - en and the interpolation 
coefficients [a^, . . ., a^^.^^} are obtained by solving the system 
of A^int + 1 equations: 



+ cr 



"ml 

I l<j<k<NiM 
Mm \^ 

-^Q-j for/ = l,...,A^i„ 

J=l ) 

'Win, ^^ 



cr - aQ + cr 

\j=l 1 l<j<k<Ni„, 

The system can be further simplified to: 



Mm 

Z(^ 

M 




fij) 






fMim 


-I- 


1 - 


Z 



(26) 



l<j<k<Ni„, 



where the first Mnt equations form a linear system which must 
be solved to obtain the Q'j=i,...,Ari,„ and where substituting these 
values in the last equation yields the value of q-q- It is worth 
noting that by using the covariances instead of the structure 
function, the system in Eq. (26) is equivalent to: 



Mim 

Ci,j aj - Co,; for i - I, 



M„, 



(27) 



2 2 



The expressions for the interpolations coefficients for the 
different cases illustrated by Fig. 4 are derived in Appendix A. 
To assess the accuracy of the statistics approximated by the 
fractal operator, we have computed the structure function of 
phase screens w computed by our implementation of the mid- 
point algorithm, i.e. as w - K m with u ~ MO,T). Fig- 
ure 5 shows that the 2D structure function is almost isotropic 
and demonstrates good agreement of our approximation to the 
theoretical law. 

D. The inverse operator 

According to the factorization in Eq. (16), the inverse of K is: 
= K„' ■ . . . ■ ' ■ K7' . (28) 




-200 2O0 100 200 300 

abscissae (units: r ) distance (units: /- ) 



Fig. 5. Structure function. Left: 2D isocontours; right: 
ID profile computed by radial averaging. Solid lines: Kol- 
mogorov law 6.88 x (r/ro)^^^; dotted lines: average of 1000 
structure functions generated with the mid-point method. 



In section 3 B, the inverse of the outermost operator Kp has 
been derived and shown to be sparse — see Eq. (22). To com- 
pute the Ky''s for the inner scales {j < p), it is sufficient to 
solve Eq. (23) for mq, which trivially yields: 



Wo • 



Y 



(29) 



where {w\ , . . . , WNi^^} are the neighbors of wo (Fig. 4). Since in 
Eq. (29), the My's only depend on the w/s, the Kt''s can be ap- 
plied in any order. However, by proceeding from the smallest 
scales toward the largest ones as in Eq. (28), the operator K"' 
can be performed in-place. This property may be important to 
avoid memory page faults and to speed-up the computation. 
Finally, from Eq. (21) and Eq. (29), it is clear that applying 
the Kj''s requires exactly as many operations as for the K/s 

and that computing K ' ■ u requires 0{N) operations. 

E. The transpose operator 

Iterating from the smallest scale to the largest one, it is easy 
to derive an algorithm to apply the transpose operator = 
Kj^ ■ . . . ■ Kt ■ K| to a given vector. The following algorithm 
computes z - ■ v for any input vector v: 



copy input vector: z <— v 
from the smallest scale to the largest scale, do 
fory = 1, . . .,Mnt do 

done 

zo <- Q'O Zq 
done 

apply Kjim at the largest scale of z 
return z 

It is important to note that the loop must be performed in- 
place for the algorithm to work. From the structure of this 
algorithm, it is clear that the multiplication of a vector by the 
transpose operator is performed in 0{N) operations. 

F. The inverse transpose operator 

The operator - Kj^^ ■ ■ . . . ■ K^^ works from the largest 
scale to the smallest one. The following algorithm computes 
z = K"^ ■ V for any input vector v: 
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copy input vector: z <— v 

apply Kj^^j at the largest scale of z 

from the largest scale to the smallest scale, do 

zo <- ^o/ao 

for; = l,...,Mntdo 

done 
done 
return z 

Again, the operation can be done in-place (the copy of the in- 
put vector V is only required to preserve its contents if needed), 
and the number of operations is 0(N). 

4. Preconditioning 

Preconditioning is a general means to speed up the conver- 
gence of iterative optimization methods [ 1 5] such as the PCG 
algorithm described in Fig. 1. Preconditioning is generally in- 
troduced as finding an invertible matrix M such that the spec- 
tral properties of M"' ■ A are more favorable than that of A 
(i.e. lower condition number and/or more clustered eigenval- 
ues), and then the transformed system 



A x ^M^ b 



(30) 



which has the same solution as the original system A ■ x - b 
can be solved in much fewer iterations. In this section, 
we consider different means for preconditioning the phase 
restoration problem: explicit change of variables and diago- 
nal preconditioners. 

A. Fractal operator as a preconditioner 

Preconditioning is also equivalent to an implicit linear change 
of variables [12]: using the preconditioner M - ■ C in 
the algorithm of Fig. 1 is the same as using the (unprecondi- 
tioned) conjugate gradient algorithm to solve the optimization 
problem with respect to x - C ■ x. Following this we have 
considered using our statistically independent modes to solve 
the problem with respect to the variables u - K ' ■ w. In 
this case, it is however advantageous in terms of the number 
of floating points operations to use an explicit change of vari- 
ables and to directly solve the problem for u rather than for w 
with a preconditioner M = K"^ ■ K"'. Introducing this change 
of variable in Eq. (7) and using Eq. (15), the system to solve 
becomes: 

(K'r ■ ■ C„' ■ S ■ K + l) ■ M = (k^ ■ S'^ ■ C„' ■ d) . (31) 

After u is found by the iterative algorithm, the restored wave- 
front is given w - K m. We expect improvements in the 
convergence of the iterative method by using u instead of w 
because this yields an a priori covariance matrix equal to the 
identity matrix [32]. Improved speedup may be still possible 
by using a preconditioner on u as we discuss in the following. 

B. Diagonal preconditioners 

Diagonal preconditioners may not be the most efficient ones 
but are very cheap to use [15] and are thus considered here. 



When the variable x in Eq. (8) follows known statistics, an 
optimal preconditioner M can be computed so that M"' ■ A is, 
on average, as close as possible to the identity matrix. This 
closeness can be measured in two different spaces: in the data 
space or in the parameter space. 

In the data space, this criterion is written: 

M = argmin<||A -x-M- x\\^) 

M 

= -2{m- A) ■ {x ■ X ) 



dm 

M C;, = A C;, , 



(32) 



where Cj: = {x ■ x^) is the covariance matrix of jc. Of course, 
if M is allowed to be any matrix and since has full rank, the 
solution to Eq. (32) is M = A. However, for a diagonal pre- 
conditioner, M = diag(ni), only the diagonal terms of Eq. (32) 
have to be considered; this yields: 



M = diag(m) = diag(A ■ C;,) ■ diag(C;,)- 



(33) 



For X - u ~ A/(0, 1) then = I and Eq. (32) simplifies to: 

M = diag(A) , (34) 

which is the well known Jacobi preconditioner [15]. 

Taking Q = M"' and minimizing the statistical distance in 
the parameter space yields: 



Q - arg min(||Q ■ A - x - x 
Q 

5<||Q-A-^-^||2) 







Q A 



2(Q-A-I)-C;,-A' 



(35) 



For a diagonal preconditioner, Q - diag(5), only the diagonal 
terms of Eq. (35) have to be considered; hence: 

Q = diag(9) = diag(C, ■ AT) ■ diag(A ■ C, ■ A^)-' . (36) 

Finally, when x - u ~ MO, I): 



Qu 



and 



0. 



(37) 



In contrast to the Jacobi preconditioner, the optimal precondi- 
tioner Q is expensive to compute since every element of ma- 
trix A must be evaluated to evaluate the denominator. This 
however has to be done only once and for all for a given left- 
hand side matrix A. The improvements given by the diagonal 
preconditioners in Eq. (34) and Eq. (37) are compared in the 
next section. 

5. Simulations and Results 

A. Summary of the various possibilities 

Our previous study gives rise to 6 diflferent possibilities to 
solve Eq. (8). The first method is based on Eq. (15) to iter- 
atively solve for w. 



(S^ ■ C„' ■ S + . K-') ■ H- = ■ C„ 



(38) 
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using the sparse model matrix S and the fractal operators 
and K"^ introduced in Sect. 2E and Sect. 3. Although the a 
priori covariance matrix of w is not the identity, we have tried 
two other methods by assessing the speedup brought by each 
of the two diagonal preconditioners defined in Eq. (34) and 
Eq. (37), with A = ■ C„' ■ S + R-t . K"'. 

Solving the problem in our statistically independent modes 
corresponds to a forth method, requiring to iteratively solve; 

(k^ ■ S'^ ■ C„' ■ S ■ K + 1) ■ H = ■ ■ C„' ■ (39) 

for u and then do w = K m. For the two last methods, we 
use with Eq. (39), one of the two preconditioners defined in 
Eq. (34) and Eq. (37) with A = ■ ■ €„' ■ S ■ K + 1. In this 
case, C„ = I so we expect somewhat faster convergence. 

B. Comparison of the rates of convergence 

When comparing the efficiency of the six different possibili- 
ties, we need to take into account that the number of floating 
point operations may be different for each of them. The aim is 
not to derive an accurate number of operations which would 
depend on the specific implementation of the algorithms, but 
rather to get a general estimate. For instance, the dependence 
of the K on ro can be factorized out and included in operator 
C„ with no extra computational cost. This kind of optimiza- 
tion was not considered here. As detailed in Appendix B, the 
number of operations is marginally increased by the precon- 
ditioning and does not depend on which variables {w or u) are 
used when starting from an arbitrary initial vector A small 
difference only appears when starting the algorithms with an 
initial zero vector, as summarized in Table 1 . 

For wavefront reconstruction, when comparing the total 
number of operations, A^ops, for a given number of (P)CG iter- 
ations, Mtei-, such that A^iter > 1, we will use these equations: 

A^CG ~ (^overhead + 33 A^^'g) N , 

(40) 

A^PCG ~ (A^overhead + 34 A^^'^'q) A? , 

where A^overhead = 4 whcn working with variable w, and 
A^oveihead = 10 whcn explicitly working with variable u. 

In order to assess the speed of the reconstruction, we have 
tested the different wavefront reconstruction algorithms on 
a number of different conditions. For every simulation, the 
wavefront sensor sampling is such that the size of the Shack- 
Hartmann subaperture is equal to Fried parameter ro. A wave- 
front is first generated by applying the fractal operator K to 
a vector of normally distributed random values like in sec- 
tion 3 C. The measurements are then estimated using the cur- 
rent wavefront sensor model, S, and a stationary uncorrected 
random noise n is added to the simulated slopes in accordance 
with Eq. (1). The noise level is given by its standard deviation 
Cnoise in radians per subaperture, where the radians here corre- 
spond to phase differences between the edges of the subaper- 
tures. At each iteration of the algorithm, the residual wave- 
front is computed as the difference between the current solu- 
tion and the initial wavefront. The root mean squared error 
of the residual wavefront is computed over the pupil, piston 




10 10 10 10 



number of floating point operations 

Fig. 6. Phase error as a function of the number of opera- 
tions. Curves are the median value of 100 simulations with 
Z)/ro = 65, (Tnoise = 1 rad/subapcrturc and ro has the same size 
as one subaperture. Solid curves are for CG, dashed curves 
are for PCG with Jacobi preconditioned dotted curves are for 
PCG with optimal diagonal preconditioner Thin curves are 
for (P)CG onto the wavefront samples w, whereas thick curves 
are for (P)CG onto the wavefront generator u. 

removed. The piston mode is the only removed mode. A cen- 
tral obscuration is always introduced, with a diameter 1/3 the 
diameter of the pupil. 

The graphs presented are for two AO system of size 65 x 65 
(c/ Figures 6, 7, 8, and 9) and 257 x 257 (cf Figures 10, 
and 11). Several levels of noise from 1 rad/subaperture down 
to 0.05 rad/subaperture are examined. They correspond to the 
levels of photon noise obtained with ~ 7 to ~ 3000 detected 
photons per subaperture. On each curve, the 6 algorithms are 
compared. All the curves plot the median value obtained for 
100 simulations under the same conditions. The different al- 
gorithms were applied to the same simulated wavefronts and 
sensor data. The graphs have been plotted assuming a number 
of floating point operations given by Eq. (40), where here the 
number of unknowns isN = 4225 and A^ = 66049 for AO sys- 
tems 65 X 65 and 257 x 257 respectively. Various observations 
can be drawn from these curves as dicussed in what follows. 

Solving by using w as unknowns is much slower than using 
H, by more than one order of magnitude for a 65 x 65 system, 
and 2 orders of magnitude for 257 x 257. This demonstrates a 
stunning efficiency for the fractal operator used as a precondi- 
tioner. With w, the algorithm does not show any improvement 
of the residual error for a long time before finding its way 
toward the solution. In contrast, the very first steps with u al- 
ready show a tremendous reduction of the residual error. For 
instance, this feature is critical if the number of iterations is to 
be limited to a fixed value as could be the case in closed-loop. 

Using Jacobi or optimal diagonal preconditioners has not 
the same effect when working in h" or in m space. When solv- 
ing for w, the preconditioners are only useful at the very end of 
the convergence, mainly in the case of high signal-to-noise ra- 
tio. Thus they are not very helpful to reduce the computational 
load. In contrast, the effect of the diagonal preconditioners is 
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number of floating point operations 

Fig. 7. Same as Fig. 6 but for cr„ai^f. - 0.5 rad/subaperture. 
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number of floating point operations 

Fig. 8. Same as Fig. 6 but for cTnoise =0.1 rad/subaperture. 

very effective from the beginning when working with u. We 
may notice that the difference between the two diagonal pre- 
conditioners is significant but not critical. The optimal diago- 
nal preconditioner yields slightly faster convergence. 

When CTnoise decreases, the convergence of the two fastest 
methods takes longer to reach a lower level of residual errors 
but the rate of convergence keeps steady. This is analyzed in 
more detail in the next section. 

C. Number of iterations 

From the previous section, we now consider only the fastest 
method, using both u as unknowns and the optimal diagonal 
preconditioner. The aim here is to assess the number of it- 
erations needed to restore the wavefront. As in the previous 
section, we consider one subaperture per ro, so the variance of 
the incoming wavefronts increases with the size of the system. 
Figure 1 2 shows how the residual phase variance decreases at 
each iteration for various configurations of the system, in size 
(33 X 33, 65 X 65, 129 x 129, 257 x 257), and in noise level 
(cr- . - 1, 0.09 and 0.01 rad^/rn). In the first iterations, we 

^ noise / uy 

can see that the behavior of the algorithm does not depend on 
the signal to noise ratio. In contrast, the final value obtained 




number of floating point operations 



Fig. 9. Same as Fig. 6 but for cr^oise - 0.05 rad/subaperture. 



100.0 



1.0 




10*** 10*' 10*" 10*' 10*" 



number of floating point operations 

Fig. 10. Same as Fig. 6 but for D/ro - 257 and o-noise - 1 
rad/subaperture. 



100.0 




10*° 10*' 10*' 10*" 10*'° 

number of floating point operations 

Fig. 11. Same as Fig. 6 but for D/ro = 257 and o-„oise = 0.5 
rad/subaperture. 
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Fig. 12. Decrease of the residual phase variance as a function 
of the number of iterations when using u as unknowns and op- 
timal diagonal preconditioner Each curve is the median value 
of 100 simulations. Three sets of curves are plotted for differ- 
ent values of cr^^.^^- 1 (soUd), 0.09 (dashed), and 0.01 rad^/ro 
(dotted). In each set of curves, the size of the system increases 
from bottom to top: 32, 64, 128 and 256 subapertures along 
the diameter of the pupil. Levels of Strehl ratios are indicated. 
The curves show that 5 to 10 iterations are enough in most 
cases for a full reconstruction. 
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number of iterations 



Fig. 13. The same curves as those in Fig. 12 are plotted here, 
normalized by the initial variance of the phase. This shows 
a high relative attenuation (~ 1/40) after the first iteration, 
in any configuration. In each set of curves: cr^^^^^ - 1 (solid), 
0.09 (dashed), and 0.01 rad'/ro (dotted); the size of the system 
increases from top to bottom: 32, 64, 128, and 256 subaper- 
tures along the diameter of the pupil. 

does not depend on the size of the system. Strehl levels cor- 
responding to the residual phase variance are indicated. The 
curves show that, whatever the size of the system, only 5 to 10 
iterations are enough for a reconstruction starting from zero. 

In order to remove the effect of starting from different ini- 
tial phase variances, the same curves have been normalized 
by the initial variance on Fig. 13. We can see that the de- 
scent of FRiM follows the same path for all the simulations. 



and is stopped at different values of the final variance, which 
depends on the signal to noise ratio. Along this path, the vari- 
ance is already reduced by a factor ~ 1 /50 at the first iteration, 
~ 1/170 at the second iteration and more than ~ lO"'* at iter- 
ation 6. This steep descent will be an asset in closed-loop. 

6. Conclusion 

We have introduced FRiM, a new minimum variance iterative 
algorithm for fast wavefront reconstruction and fast control of 
an adaptive optics system. Combining fast regularization and 
efficient preconditioning, regularized wavefront reconstruc- 
tion by FRiM is an 0{N) process, where is the number of 
wavefront samples. 

FRiM takes advantage of the sparsity of the model matrix S 
of wavefront sensors (or interaction matrices) and makes use 
of a "fractal operator" K for fast computation of the priors. 
Based on a generalization of the mid-point algorithm [23], 
K is not sparse but is implemented so that it requires only 
0{N) ^ 6N operations. Our modifications with respect to the 
original algorithm allow the operator to be invertible and the 
generated wavefront to be stationary. We have derived algo- 
rithms for computing K and K"^ in the same number 
of floating point operations. In our simulations, we consider 
a modified Kolmogorov law but any stationary structure func- 
tion or covariance can be implemented in our approach. The 
property of stationarity is expected to be helpful for turbulence 
tomography. 

Another breakthrough comes from the efficiency of the 
fractal operator when used as a preconditioner Combining 
a fractal change of variables and an optimal diagonal precon- 
ditioner, we were able to reduce the number of iterations in the 
range of 5 - 10 for a full wavefront reconstruction whatever 
is the size of the AO system. The exact number of iterations 
mainly depends on the signal to noise ratio of the measure- 
ments. 

It is beyond this work to compare with all the other meth- 
ods currently studied in response to the huge increasing of the 
number of degrees of freedom for the AO system on ELTs. 
Nevertheless, we can easily compare to standard vector ma- 
trix multiplication (VMM). Assuming uncorrected noise, the 
simulations show that the number of operations with FRiM is 
A^ops ~ (23 + 34 A^iter) A^, where the number of PCG iterations 
is A^iter ^ 10 for any number of degrees of freedom A^. For 
up to A^ = 1.3 X 10"* de grees of freedom (i.e. D/ro < 128), 
one wavefront estimation (from scratch) involves < 6 x 10^ 
operations, that is a bandwidth of ~ 500 Hz for a machine 
capable of 3 Gflops which is typical of current workstations. 
Conversely, conventional (non-sparse) matrix multiplication 
would require ~ 4A^^ ~ 7 x 10^ operations to compute the 
wavefront: our method is more than 100 times faster. Further- 
more, since the operations can be done in-place, it is expected 
that the computation with FRiM could all be done in cache 
memory. 

For simulating very large AO systems {e.g. atmospheric 
tomography on ELT's), the speed of the current version of 
FRiM is already an asset. For real-time control of AO sys- 
tems, FRiM algorithm can be parallelized to run on several 
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CPU's. Being an iterative method (unlike Fourier methods), 
FRiM could be used to improve the estimation of the wave- 
front from any pieces of new data as soon as it becomes avail- 
able. Hence, FRiM does not need all the measurements in a 
closed-loop system. A fast iterative method that gives inter- 
mediate results with only a part of the measurements opens the 
way to new control approaches for reducing the effect of the 
delay. A further advantage of FRiM is that it accounts for the 
statistics of the turbulence which not only yields a better esti- 
mation of the residual phase [33] but also helps to disentangle 
ambiguities such as unseen modes in atmospheric tomogra- 
phy. In this paper, we assume that the structure function is 
perfectly known. Bechet [34] has shown that it is sufficient to 
not overestimate ro by more than a factor ~ 2 to benefit from 
the advantages of taking into the priors. 

The next step of this work is to extend the theory to closed- 
loop and to assess the performances and the properties of the 
algorithm in this regime. Since the wavefront is not allowed 
to change a lot from one step of the AO loop to the other, the 
algorithm will always starts close to the solution: the number 
of iterations is expected to be yet lower This study is not yet 
completed but preliminary results have proved the efficiency 
of FRiM in the case of closed-loop adaptive optics [35, 36]. 
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Appendix A: Derivation of the interpolation coefHcients 

In this appendix, we detail the computation of the interpola- 
tion coefficients involved in the different configurations shown 
by Fig. 4. Denoting r the step size in the grid before the re- 
finement, the distances between the points considered in this 
refinement step are: V2r, r, r/V2, or r/2 (Fig. 4). Hence the 
only covariances required in our computations are: 



(Al) 



where c{r) and f{r) are respectively the covariance and the 
structure function for a separation r. 

1. Square configuration 

For the interpolation stage illustrated by the top-left part of 
Fig. 4 and according to Eq. (27), the interpolation coefficients 



{o-i, Q'2, a3, q;4) are obtained by solving: 
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Solving this linear system and plugging the solution into 
Eq. (27) leads to: 



a\ — a2 — ctT, — Q'4 



C2 



Co H- 2 C3 -I- C4 ' 
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Co H- 2 C3 -H C4 



Note that the sign of ao is irrelevant. 
2. Triangle configuration 

In original mid-point algorithm [23], the values at the edges 
of the support (top-right part of Fig. 4) were generated from 
only the two neighbors on the edge, ignoring the third in- 
terior neighbor (denoted in the figure). Here, according 
to Eq. (27), the interpolation coefficients {ai,a2,cij,] for this 
stage are obtained by solving: 
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Solving this linear system and plugging the solution into 
Eq. (27) leads to: 



ffl = Q'2 



a'3 = 



Cl (co - C2) 
Co (co -H C3) - 2 c\ 
Cl (co - 2 C2 -H C3) 
Co (co -H C3) - 2 c] 



(A3) 



ao - 



c\ (3 Co - 4 C2 -H C3) 



Co - 



Co (co H- C3) - 2 c\ 



3. Diamond configuration 



The interpolation coefficients for the stage in the bottom part 
of Fig. 4 can be deduced from Eq. (A2) by replacing r by r/V2, 
then: 



Cl = a2 = a3 = a4 = 



ao - 



Cl 



Co + 2 C2 H- C3 



CO 



4c2 



(A4) 



Co H- 2 C2 -H C3 

Appendix B: Computational Burden 

In order to estimate the number of floating point operations, 
we need to carefully detail the steps of the CG method and 
count the number of operations involved in the multiplication 
by the different linear operators S, K, etc. . Figure 1 summa- 
rizes the steps of the (PCG) algorithm [15] to solve Eq. (8). 
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algorithm step floating point operations 



initialization: general case 


~ 25 A' 


zero initial vector in u space 


~ 12N 


zero initial vector in vc space 


~ 6N 


1st CG iteration 


~ 31N 


any subsequent CG iteration 


~ 33 Af 


total after N^^y > 1 iterations 


~(23 + 33A'i,e,-)A' 


1st PCG iteration 


~ 32N 


any subsequent PCG iteration 


~ 34 


total after A'jie, > 1 iterations 


~(23 + 34A'ue,-)A' 



Table 1 . Number of operations involved in conjugate gradi- 
ents (CG) and preconditioned conjugate gradients (PCG) ap- 
plied to the wavefront restoration problem solved by our algo- 
rithm. The integers and A^itei are respectively the number of 
unknowns and number of iterations. For a reconstruction, we 
assume an initial null guess in the initialization step: in this 
case the number of operations at this step is reduced down to 
~ 6N 01 ~ 12 N when respectively w or u are used as un- 
knowns. 

If the unknowns are the wavefront samples, then x -w and: 
A = ■ C„' ■ S + ■ , 

Starting the algorithm with Wq, the initial residuals write: 

ro = ft - A ■ H-o 

= ■ C„' ■ (rf - S ■ wq) - . ■ wq . (Bl) 

If the unknowns are the wavefront generators, then x = u 
and: 

A = KT.S'^-C„' - S-K + I, 
ft = KTS^C„'rf, 

where I is the identity matrix. Starting the algorithm with hq, 
the initial residuals are: 

ro = ft - A ■ Mo 

= ■ S'^ ■ C„' ■ (rf - S ■ K ■ Ho) - Mo . (B2) 

Making use of possible factorizations (some of the a/'s have 
the same values), applying any one of the operators K, K^, 
K or K"^ involves the same number of floating point oper- 
ations: 

A?ops(K) = A?ops(K'^) = A^op,s(K-') = A?ops(K-T) 
= 6Nu - 14 
~ 6A^, 

where A^ is the number of degrees of freedom of the system, 
Nu ~ N is the number of elements in vector u and, in our 
notation, A^ops(L) is the number of floating point operations 
required to apply a linear operator L to a vector 



Since we consider uncorrected data noise, C„' is diagonal 
and: 

A^ops (C„') = M~2A?; 

however note that these ~ 2N floating point operations per 
iteration can be saved for stationary noise (i.e. C„' oc I). 

For Fried model of wavefront sensor and after proper fac- 
torization: 

A^ops(S)=A?ops(S^)~4A?. 

This assumes, in particular, that the data were pre-multiplied 
by 2 (see Eq. (11)). 

Finally, whatever the unknown are {w or «), the total num- 
ber of floating point operations required to apply the left hand 
side matrix A to a given vector is: 

A^ops(A) ~ INopAK) + 2 A?ops(S) + A?ops(C„') + A^ 
~ 23N. 

The last A^ comes from the addition of likelihood and regular- 
ization terms. 

From equations (Bl) and (B2), using either w or « as the 
unknowns, initialization of the CG, i.e. computation of the 
initial residuals ro, involves 

A^ops(ro) ~ 2 A?ops(K) + 2 N,^,{S) + A?ops(C„>) + M + N 
~ 25 A^ 

operations. Note that, if the algorithm is initialized with 
jco = (a vector of zeroes), this number of operations is signif- 
icantly reduced down to ~ 6N and ~ 12N when respectively 
w and u are used as unknowns. Also note that there may be 
additional ~ 6N operations to compute w from u when nec- 
essary. 

Whatever are the considered variables, the number of un- 
knowns is ~ A^, hence any dot product in the CG algorithm 
involves 2N - I ~ 2N floating point operations. The first 
CG iteration (Fig. 1 ) requires two dot products (2 A^ - 1 ~ 2 A^ 
floating point operations each) to compute and ff^, apply- 
ing A once and two vector updates (involving ~ 2N opera- 
tions each); hence a total of ~ 31 A^ operations. Any subse- 
quent iteration requires an additional vector update to com- 
pute the conjugate gradient direction; hence ~ 33 A^ opera- 
tions. Finally, preconditioning by a diagonal preconditioner 
simply adds ~ A^ operations per iteration. 

The number of floating operations required by the diff'erent 
versions of the reconstruction algorithm are summarized in 
table 1 and by Eq. (40). Note that in the general case, the 
number of operations does not depend on which variables w 
or u are used. There is a difference of ~ 6N operations in 
the initialization step only when the algorithm is started with 
a zero initial vector (see table 1). 
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